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Figure 1. Distribution and two tomograms. 

1. Introduction 

The primary goal of tomography is to determine the 
internal structure of an object without cutting it, 
namely using data obtained by methods that leave the 
object under investigation undamaged. These data can 
be obtained by exploiting the interaction between the 
object and various kinds of probes including X-rays, 
electrons, and many others. After its interaction with 
the object under investigation, the probe is detected 
to produce what we call a projected distribution or 
tomogram, see Fig. 

Tomography is a rapidly evolving field for its 
broad impact on issues of fundamental nature and for 
its important applications such as the development of 
diagnostic tools relevant to disparate fields, such as 
engineering, biomedical and archaeometry. Moreover, 
tomography can be a powerful tool for many 
reconstruction problems coming from many areas of 
research, such as imaging, quantum information and 
computation, cryptography, lithography, metrology 
and many others, see Fig. 

From the mathematical point of view the 
reconstruction problem can be formulated as follows: 
one wants to recover an unknown function through 
the knowledge of an appropriate family of integral 
transforms. It was proved by J. Radon [T] that a 
smooth function f(x,y) on can be determined 
explicitly by means of its integrals over the lines in 

Let TZf{X,9) denote the integral of / along the 
line X cos 0 + y sin 9 = X (tomogram). Then 

n27r 

f{x,y) = i-A)^ J nfixcos9+ ysm9,9) (1) 

where A = is the Laplacian on and 



Figure 2. Reconstruction problems from diverse fields may be 
united within the framework of tomography. 

its square root is defined by Fourier transform (see 
Theorem [^. We now observe that the formula above 
has built in a remarkable duality: first one integrates 
over the set of points in a line, then one integrates over 
the set of lines passing through a given point. This 
formula can be extended to the iV-dimensional case by 
computing the integrals of the function / on all possible 
hyperplanes. This suggests to consider the transform 
/ I—>■ TZf defined as follows. If / is a function on 
then TZf is the function defined on the space of all 
possible {N — I)-dimensional planes in such that, 
given a hyperplane A, the value of TZf{X) is given by 
the integral of / along A. The function TZf is called 
Radon transform of /. 

There exist several important generalizations of 
the Radon transform by John [5], Gel’fand |3], 
Helgason [3] and Strichartz [S]. More recent analysis 
has been boosted by Margarita and Volodya Man’ko 
and has focused on symplectic transforms [6], on the 
deep relationship with classical systems and classical 
dynamics [ziiHiin], on the formalism of star product 
quantization [iniiiiiiii], and on the study of marginals 
along curves that are not straight lines Haul]. 

In quantum mechanics the Radon transform of 
the Wigner function [laiiiin] was considered in the 
tomographic approach to the study of quantum states 
na na and experimentally realized with different 
particles and in diverse situations. For a review 
on the modern mathematical aspects of classical and 
quantum tomography see |20j . Good reviews on recent 
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tomographic applications can be found in [2T] and 
in [22] , where particular emphasis is given on maximum 
likelihood methods, that enable one to extract the 
maximum reliable information from the available data 
can be found. 

As explained above, from the mathematical point 
of view, the internal structure of the object is 
described by an unknown function / (density), that is 
connected via an operator to some measured quantity 
g (tomograms). The tomographic reconstruction 
problem can be stated as follows: for given data 
g, the task is to find / from the operator equation 
TZf = g. There are many problems related to the 
implementation of effective tomographic techniques 
due to the instability of the reconstruction process. 
There are two principal reasons of this instability. 
The first one is the ill-posedness of the reconstruction 
problem: in order to obtain a satisfactory estimate 
of the unknown function it is necessary an extremely 
precise knowledge of its tomograms, which is in general 
physically unattainable [23]. The second reason is the 
discrete and possibly imperfect nature of data that 
allows to obtain only an approximation of the unknown 
function. 

The first question is whether a partial information 
still determines the function uniquely. A negative 
answer is given by a theorem of Smith, Solomon 
and Wagner [24], that states: “A function / with 
compact support in the plane is uniquely determined 
by any infinite set, but by no finite set of its 
tomograms”. Therefore, it is clear that one has to 
abandon the request of uniqueness in the applications 
of tomography. Thus, due to the ill-posedness of 
reconstruction problem and to the loss of uniqueness 
in the inversion process, a regularization method has 
to be introduced to stabilize the inversion. 

A powerful approach is the introduction of 
a Mumford-Shah (MS) functional, first introduced 
in a different context for image denoising and 
segmentation [25] . The main motivation is that, in 
many practical applications, one is not only interested 
in the reconstruction of the density distribution /, 
but also in the extraction of some specific features or 
patterns of the image. An example is the problem of 
the determination of the boundaries of inner organs. 
By minimizing the MS functional, one can find not 
only (an approximation of) the function but also its 
sharp contours. Very recently a MS functional for 
applications to tomography has been introduced in the 
literature [2511271 HHj. Some preliminary results in this 
context are already available but there are also many 
interesting open problems and promising results in this 
direction, as we will try to explain in the second part 
of this article. 

The article is organized as follows. Section [^ 
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contains a short introduction to the Radon transform, 
its dual map and the inversion formula. Section [^ is 
devoted to a brief discussion on the ill-posedness of the 
tomographic reconstruction and to the introduction of 
regularization methods. In Section [^ a MS functional 
is applied to tomography as a regularization method. 
In particular, in Subsection |4.1| the piecewise constant 
model and known results are discussed together 
with a short list of some interesting open problems. 
Finally, in Section [^ we present an electrostatic 
interpretation of the regularization method based on 
the MS functional, which motivates us to introduce an 
improved regularization method, based on the Blake- 
Zisserman functional [^, as a relaxed version of the 
previous one. 

2. The Radon transform: definition and 
inversion formula 

Consider a body in the plane M ?, and consider a beam 
of particles (neutrons, electrons. X-rays, etc.) emitted 
by a source. Assume that the initial intensity of the 
beam is Jq. When the particles pass through the body 
they are absorbed or scattered and the intensity of the 
beam traversing a length As decreases by an amount 
proportional to the density of the body /r, namely 

AI/1 ——As fi{s), (2) 

so that 

I{s) = Iq exp J p,{r) dr^ . (3) 

A detector placed at the exit of the body measures the 
final intensity J(s) and then from 

— In = f /i(r) dr (4) 

4 Jo 

one can record the value of the density fi integrated 
on a line. If another ray with a different direction is 
considered, with the same procedure one obtains the 
value of the integral of the density on that line. 

The mathematical model of the above setup is the 
following: Given a smooth function f{x) on the plane, 
a; S and a line A, consider its tomogram, given by 

dm{x), (5) 

where dm is the Euclidean measure on the line A. In 
this way, we have defined an operator TZ that maps a 
smooth function / on the plane into a function TZf 
on P^, the manifold of the lines in 

We ask the following question: If we know the 
family of tomograms {'R-f{X))x^p 2 , can we reconstruct 
the density function fl The answer is affirmative and 
in the following we will see how to obtain this result. 

Let us generalize the above definitions to the case 
of an A-dimensional space. Let / be a function defined 
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Figure 3. Parametrization of the hyperplane A using its signed 
distance X from the origin and a unit vector ^ perpendicular to 

A. 

on integrable on each hyperplane in and let 
be the manifold of all hyperplanes in . The 
Radon transform of / is defined by Eq. ([^, where dm 
is the Euclidean measure on the hyperplane A. Thus we 
have an operator TZ, the Radon transform, that maps 
a function / on into a function TZf on P^, namely 
/ I—> TZf. Its dual transform, also called back projection 
operator, g i—>■ Tg associates to a function g on P^ the 
function Ig on given by 

= [ gW d/r(A), (6) 

J xGX 

where d^ is the unique probability measure on the 
compact set {A S P^|a; G A} which is invariant under 
the group of rotations around x. 

Let us consider the following covering of P^ 

K X ^ P^, (X, 0 A, (7) 

where is the unit sphere in Thus, the 

equation of the hyperplane A is 

A = {xeR^IX-^-x = 0}, (8) 

with a ■ b denoting the Euclidean inner product of 
a,b G . See Fig. Observe that the pairs 

(X, ^), (—X, —^) G M X are mapped into the 

same hyperplane A G P^. Therefore Q is a 
double covering of P^. Thus P^ has a canonical 
manifold structure with respect to which this covering 
mapping is differentiable. We identify continuous 
(differentiable) functions g on P^ with continuous 
(differentiable) functions g on M x satisfying 

g{x,0 = g{.-x,-0. 

We will momentarily work in the Schwartz space 
S{R^) of complex-valued rapidly decreasing functions 



Figure 4. X(7Zf){x) is the potential at x generated by the 
charge distribution /. 


on In analogy with S{R^) we define 5(M x S^“^) 
as the space of C°° functions g on K x which for 
any integers m > 0, any multiindex a G N^, and any 
differential operator U on satisfy 


sup 

AfGR.fGS™-! 




< -1-00. 


(9) 


The space 5(P^) is then defined as the set of g G 
5(M X S^"i) satisfying g^-X, -^) = g{X, ^). 

Now we want to obtain an inversion formula, 
namely we want to prove that one can recover a 
function / on from the knowledge of its Radon 
transform. In order to get this result we need a 
preliminary lemma, whose proof can be found in |20] . 
which suggests an interesting physical interpretation. 
Lemma 1. Let f G S(R-^) andV(x) = l/|a;|, x G R^, 
X ^ 0 . Then 


T{TZf){x) = aN f *V, 


( 10 ) 


where qm depends only on the dimension N, and * 
denotes the convolution product, 

{f*g)(,x)=[ f{y)gix-y) dy. (11) 

A physical interpretation of Lemma is the 
following: if / is a charge distribution, then the 
potential at the point x generated by that charge is 
exactly X(77./(a;)), see Fig.[^ Notice, however, that the 
potential of a point charge scales always as the inverse 
distance independently of the dimension N, and thus 
it is Coulomb only for N = 3. The only dependence 
on N is in the strength of the elementary charge qn- 
This fact is crucial: indeed, the associated Poisson 
equation involves an A^-dependent (fractional) power of 
the Laplacian, which appears in the inversion formula 
for the Radon transform. 


Theorem 1. Let f G 5(]R^). Then 

where (—A)“, with a > 0, is a pseudodifferential 
operator whose action is 


{{-Arf){x) 


|fc|2“/(A:)e‘'='^ 


dfc 

(2^’ 


(13) 
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where f is the Fourier transform of f, 
f{k)= ( f{x)e-''^--Ax. 


(14) 


The proof of Theorem can be found in m- 
Equation (12 1 says that, modulo the final action of 
(—A)t^, the function / can be recovered from its 
Radon transform TZf by the application of the dual 
mapping T: first one integrates over the set of points 
in a hyperplane and then one integrates over the set of 
hyperplanes passing through a given point. Explicitly 
we get 

f(^)= / 7^/(^x,e)de, (15) 

2(27rj Jsw-i 

which has the following remarkable interpretation. 
Note that if one fixes a direction ^ G then 

the function TZf{£, ■ x,f,) is constant on each plane 
perpendicular to i.e. it is a (generalized) plane wave. 
Therefore, Eq. (151 gives a representation of / in terms 


of a continuous superposition of plane waves. A well- 
known analogous decomposition is given by Fourier 
transform. When N = 2, one recovers the inversion 
formula ([^ originally found by Radon [T] . 


3. Instability of the inversion formula with 
noisy data 

We have defined the Radon transform of any function 
/ G as TZf. The following theorem a contains 

the characterization of the range of the Radon linear 
operator TZ and the extension of TZ to the space of 
square integrable functions L^(K^). 

Theorem 2. The Radon transform TZ is a linear one- 
to-one mapping of S{M.^) onto Sh{V^), where the 
space SniV^) is defined as follows: g G Sh(V^) if 
and only if g G S{¥^) and for any integer fc G N the 
integral 

[g{X,OX'^dX (16) 

Jm. 

is a homogeneous polynomial of degree k m ^i,... 
Moreover, the Radon operator TZ can be extended to a 
continuous operator from L^(IR^) and L^(P^). 

In medical imaging, computerized tomography 
is a widely used technique for the determination of 
the density / of a sample from measurements of 
the attenuation of X-ray beams sent through the 
material along different angles and offsets. The 
measured data g are connected to the density / via 
the Radon Transform TZ. To compute the density 
distribution / the equation g = TZf has to be inverted. 
Unfortunately it is a well known fact that TZ is not 
continuously invertible on L^(P^) |4], and this imply 
that the problem of inversion is ill-posed. For this 


reason, regularization methods have to be introduced 
to stabilize the inversion in the presence of data noise. 

We discuss ill-posed problems only in the 
framework of linear problems in Hilbert spaces |30j . 
Let JF be Hilbert spaces and let A be a linear 
bounded operator from into JU. The problem 

given g G JU, find f G such that Af = g (17) 

is called well-posed by Hadamard (1932) if it is 
uniquely solvable for each g G JF and if the solution 
depends continuously on g. Otherwise, 0 is called 
ill-posed. This means that for an ill-posed problem the 
operator A~^ either does not exist, or is not defined on 
all of or is not continuous. The practical difficulty 
with an ill-posed problem is that even if it is solvable, 
the solution of Af = g need not be close to the solution 
of Af = g’^ if g'^ is close to g. 

In general A~^ is not a continuous operator. 
To restore continuity we introduce the notion of a 
regularization of A~^. This is a family (T.y).y>o of linear 
continuous operators T^ : —>■ which are defined 

on all GifG and for which 

hrn^ T-^g = A~'^g (18) 

on the domain of A~^. Obviously ||T.y|| —>■ -boo as 
7 —)■ 0 if A~^ is not bounded. With the help of a 
regularization we can solve 0 approximately in the 
following sense. Let g'^ G be an approximation to 
g such that \\g — g'^|| < e . Let 7 (e) be such that, as 
e —>■ 0, 

7 (e) ^0, ||T^(,)||e^0. (19) 

Then, as e —>■ 0, 

- A-I 5 II < ||r^(,)(/-ff)|! 

+ ||(r^(,)-A-i)5|| 

< ||r^(,)||e+||(r^(,)-A-i)g|| ^0. 

( 20 ) 

Hence, Ty(e) 5 '^ is close to A~^g if g"^ is close to g. 
The number 7 is called a regularization parameter. 
Determining a good regularization parameter is one of 
the crucial points in the application of regularization 
methods. 

There are several methods for constructing 
a regularization as the truncated singular value 
decomposition, the method of Tikhonov-Phillips or 
some iterative methods m- In the following section 
we present a regularization method based on the 
minimization of a Mumford-Shah type functional. 

4. Mumford-Shah functional for the 
simultaneous segmentation and reconstruction 
of a fnnction 

In many practical applications one is not only 
interested in the reconstruction of the density 
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distribution / but also in the extraction of some specific 
features within the image which represents the density 
distribution of the sample. For example, the planning 
of surgery might require the determination of the 
boundaries of inner organs like liver or lung or the 
separation of cancerous and healthy tissue. Segmenting 
a digital image means finding its homogeneous regions 
and its edges, or boundaries. Of course, the 
homogeneous regions are supposed to correspond to 
meaningful parts of objects in the real world, and 
the edges to their apparent contours. The Mumford- 
Shah variational model is one of the principal models 
of image segmentation. It defines the segmentation 
problem as a joint smoothing/edge detection problem: 
given an image g{x), one seeks simultaneously a 
“piecewise smoothed image” u{x) with a set F of 
abrupt discontinuities, the “edges” of g. The original 
Mumford-Shah functional [53], is the following: 

MS{T,u) = \\u- gWl^jy. + a f |Vu(a::)p dx 

Jd\t 

+ ( 21 ) 

where 

• D C is an open set {screen)-, 

• F C is a closed set {set of edges)-, 

• u : D ^ M. {cartoon)-, 

• Vm denotes the distributional gradient of u; 

• g G L‘^{D) is the datum {digital image)-, 

• a,/3 > 0 are parameters {tuning parameters)-, 

• denotes the {N — l)-dimensional Hausdorff 
measure. 


The squared distance in (21) plays the role 
of a fidelity term: it imposes that the cartoon u 
approximate the image g. The second term in the 
functional imposes that the cartoon u be piecewise 
smooth outside the edge set F. In other word this 
term favors sharp contours rather than zones where a 
thin layer of gray is used to pass smoothly from white 
to black or viceversa. Finally the third term in the 
functional imposes that the contour F be “small” and 
as smooth as possible. What is expected from the 
minimization of this functional is a sketchy, cartoon¬ 
like version of the given image together with its 
contours. See Fig. 

The minimization of the MS functional represents 
a compromise between accuracy and segmentation. 
The compromise depends on the tuning parameters 
a and j5 which have different roles. The parameter 
a determines how much the cartoon u can vary, if a 
is small some variations of u are allowed, while as a 
increases u tends to be a piecewise constant function. 
The parameter {3 represents a scale parameter of the 
functional and measure the amount of contours: if /3 



Figure 5. Left: image of an eye (g). Center: contours of the 
image in the Mumford-Shah model (edges F). Right: piecewise 
smooth function approximating the image (cartoon u) mi¬ 


ls small, a lot of edges are allowed and we get a fine 
segmentation. As {3 increases, the segmentation gets 
coarser. For more details on the model see the original 
paper [53], and the book [5T] . 


The minimization of the MS functional in (21) is 


performed among the admissible pairs (F, u) such that 
F is closed and u G C^{D\T). It is worth noticing 
that in this model there are two unknowns: a scalar 
function u and the set F of its discontinuities. For this 
reason this category of problems is often called “Free 
Discontinuities Problem”. Existence of minimizers of 


the MS functional in (21) was proven by De Giorgi, 
Carriero, Leaci in |33] in the framework of bounded 
variation functions without Cantor part (space SBV) 
introduced by Ambrosio and De Giorgi in |34j . Further 
regularity properties for optimal segmentation in the 
Mumford-Shah model were shown in [331 ESI EH EH]. 

Here we present a variation of the MS functional, 
adapted to the inversion problem of the Radon 
transform. More precisely, we consider a regularization 
method that quantifies the edge sets together with 
images, i.e. a procedure that gives simultaneously 
a reconstruction and a segmentation of / (assumed 
to be supported in D C directly from the 

measured tomograms g, based on the minimization of 
the Mumford-Shah type functional 


a / |V/(x)f da; -I- {3% 

Jd\t 


N-l 


(F). (22) 


The only difference between the functionals MS and 
Jms is the first term, i.e. the fidelity term, that 
ensures that the reconstruction for / is close enough 
to a solution of the equation TZf = g, whereas the 
other terms play exactly the same role explained for 
the functional MS. As explained above, in addition to 
the reconstruction of the density f, we are interested 
in the reconstruction of its singularity set F, i.e. the 
set of points where the solution / is discontinuous. The 
main difference with respect to the standard Mumford- 


Shah functional (21) is that we have to translate the 


information about the set of sharp discontinuities of g 
(and hence on the space of the Radon transform) into 
information about the strong discontinuities of /. 
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4.I. The piecewise constant model 


Here we will review the results obtained by Ramlau 
and Ring |26j concerning the minimization of (22) 


restricted to piecewise constant functions /, and then 
consider some interesting open problems. For medical 
applications, it is often a good approximation to 
restrict the reconstruction to densities / that are 
constant with respect to a partition of the body, as 
the tissues of inner organs, bones, or muscles have 
approximately constant density. 

We introduce the space PCm{D) as the space 
of piecewise constant functions that attain at most 
m different function values, where D is an open and 
bounded subset of . In other words, each / S 
PCm{D) is a linear combination of m characteristic 
functions xo.^ of sets {^k)k=i,...,m which satisfy 


m 

XOfc = XD a.e. 

fc=i 


We assume that the ftk’s are open relatively to D and 
we set Ffe = dVtk for the boundary of Tlk with respect 
to the topology relative to the open domain D. In 
this situation the edge set will be given by the union 
of the boundaries of H^’s. For technical reasons it is 
necessary to assume a nondegeneracy condition on the 
admissible partitions of D: 

Pt = {flk)k=i,...,m is admissible if > S, (23) 

for some S > 0, for all fc = 1,..., to, where denotes 
the Lebesgue measure on . 

It turns out to be convenient to split the 
information encoded in a typical function. 




(24) 


k=l 


into a “geometrical” part described by the m-tuple of 
pairwise disjoint sets = i^k)k=i,...,m which cover 
D up to a set of measure zero and a “functional” 
part given by the TO-tuple of values / = ifk)k=i,...,m- 
We also use the notation F = (Tk)k=i,...,m, for the 
boundaries F^ = dflk of flk- 

As usual when dealing with inverse problems, 
we have to assume that the data g are not exactly 
known, but that we are only given noisy measured 
tomograms g'^ of a (hypothetical) exact data set g with 
\\f - g||L2(p«) < £• 


If we restrict the functional (22) to functions in 


PCm{D) we obtain that the second term (involving 
the derivatives of /) disappears, therefore it remains 
to minimize the functional 


Jg{nj) = +/3^H^-i(Ffc), (25) 

k=l 

over PCm{D), with respect to the functional variable f 
(a vector of to components) and the geometric variable 


n (a partition of the domain D with at most to distinct 
regions satisfying the non degeneracy condition (23)). 
So the problem is to find / G PCmiD) such that 


/ — fkXn,, G PCm{D), 

k=l 

where 

(fi, /) = arg min J^(r, /). 


(26) 


(27) 


It is clear that / will depend on the regularization 
parameter (3 and on the error level e. 

Now we can state the results concerning the 


functional Jp in (25). There are several technical 
details necessary for the precise statement and proof 
of the theorems, for which we refer to the original 
paper |26j . Here we will give a simplified version of the 
theorems with the purpose of explain the main goal, 
without too many technical details. The first result is 
about the existence of minimizers of the functional Jp 
in (|2|). 

Theorem 3. For all g’^ G L^(P-^) there exists a 


minimizer {flp, fp) of the functional Jp in (25), with 
/3 > 0. 

The second result regards the stable dependence 


of the minimizers of the functional Jp in (251 on the 
error level e. 


Theorem 4. Let (<7'^")nGN be a sequence of functions 
in L^(P'^) and let g'^ G L^(P^). For all n S N, 
let {FTp , fp^) denote the minimizers of the functional 
Jp with initial data g^". If g'^" —>■ g'^ in as 

n —>■ + 00 , then there exists a subsequence of , /^") 
such that 

as j —>■ + 00 , and [FVp, f^) is a minimizer of Jp with 
initial data g'^. Moreover, the limit of each convergent 
subsequence of ($7^",/^") is a minimizer of Jp with 
initial data g^. 

Finally the last Theorem is a regularization result. 
Theorem 5. Let f* G PCmiD) be given. 


f — fkXGll: 

fc=l 

and let g* = TZf*. Assume we have noisy data 
ge g +2(pAr) ye _ < e. Let us 

choose the parameter (3 = I3{e) satisfying the conditions 
/3(e) —>■ 0 and e^/,5(e) —)■ 0 as e —>■ 0. For any 
sequence e„ —>■ 0, let in-,r) denote the minimizers 
of the functional Jp{e„) initial data g'^^ and 

regularization parameter (3 = /3(e„). Then there exists 
a convergent subsequence of iFC,f'^). Moreover, for 









Tomography: mathematical aspects and applications 


8 


every convergent subsequence with limit ($1, /) the 
function 




fc=i 


Is the nondegeneracy condition (231 necessary? 


• Can one find an a priori optimal value for the 
number m of different values? 

• Is it possible to give an a priori estimate on the 
L°°-norm of the solution (maximum principle)? 

• And finally, it would be very important for 
applications to prove the existence of minimizers 
of the functional Jms not restricted to piecewise 
constant functions /. 

We observe that all these problems are quite natural, 
and have been completely solved in the case of the 
standard Mumford-Shah functional MS in (21), see 

e.g. [211 [3S]- 

5. Electrostatic interpretation of Jms 

In this section we restrict our attention to the 
3-dimensional case. We propose an electrostatic 
interpretation of the regularization method based on 
the functional Jms discussed in the previous section. 
The intent is to give a physical explanation of the 
fidelity term ||??-/ —g||| 2 (p 3 ) in the functional that 
provide the intuition for an improved regularization 


method. For N = 3, the inversion formula (12) and the 
electrostatic identity ( |l 0 | particularize, respectively, as 
follows: for all / S 5(^) one gets 

f = (28) 


2(27r)2 


and 


X(7^/) = asf * E, (29) 

where V{x) = l/|a:| and 03 is a constant. We present 
two preliminary Lemmas. 

Lemma 2. For all real valued f G S(U^) one has 
1 


m\\ = 


2(27r)2 


|VZ(7^/)r. 


(30) 


Proof. We know that / = — 2 ( 5 ^ AI(7?./), therefore 


IXS2 


is a solution of the equation TZf = g* with a minimal 
perimeter. Moreover if f* is the unique solution of this 
equation then the whole sequence converges 

r^f\ mL\D), 

when n —>■ -l-oo. 

Finally, let us list some open problems in this 
context: 


\Tlf{X,0\^dXd^= [ f{x)I{nf){x)dx 
[-AI(7^/)](x)Z(7^/)(al) dx 


Im 2(27r)^ 
1 f 


2(2 


TT 


|VZ(7^/)(x)pdx. 


Lemma 3. For all real valued f G 5(M^) define 
1 


E = - 


and 


F = 


2(27r)- 


^VI(7^/) 




2(27r)2' 

Then 

f = X ■ E = -Aip. 


□ 

(31) 

(32) 

(33) 


Proof. 
V-E= - 


1 


• VI(7^/) = -A(^ 


2(27r) 


2(27r)2 

^ (-A)I(7^/) = /, 


where we used the inversion formula (28). 


□ 


Now we consider a measured tomogram g : ^ 
M and let us assume that g = 72,/o for some /o : t 

K. By Lemma [2]|^ it follows immediately that the 
fidelity term ||72./— 5 |j^ 2 (p 3 ) can be rewritten as follows: 

m-g\\l,^p,)=2{27rf\\E-Eg\\lHm) 


= 2(27r)2||V(/3-Vv3, 


9llL2( 


(34) 


where 
E = -- 


iVI(72/), 


1 


2(27r)2 


2(27r)2 

are the corresponding electric fields, while 


VI(g) (35) 


F = 




1 


jUg) 


(36) 


2(27r)2 ' 2(27r)2- 

are the corresponding potentials. 

With respect to the standard Mumford-Shah 


functional MS in ( 21 ), the new fidelity term in the 


functional Jms in ( 22 ) controls the distance between 


the Radon transform of / and the tomographic data 
g. The relevant difference with respect to the original 
functional is that the function / and its Radon 
transform TZf are defined in different spaces. Let 
us try to interpret the fidelity term ||72,/ — g||| 2 (p 3 ) 
from a physical point of view. A key ingredient for 
this goal is the electrostatics formulation of the Radon 
transform. This formulation can be summarized as 
follows: if we consider, in dimension 3, a function /, 
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we can think at it as a charge distribution density; if 
we apply to / first the Radon operator TZ and then its 
adjoint I we obtain, up to a constant, the electrostatic 
potential generated by the charge distribution /. This 
formulation can be stated in any dimension N: the 
difference with general potential theory in dimension 
N is that, in tomography, the potential produced by a 
point charge always scales like l/|a;|, which is the case 
of electrostatic potential only in dimension 3. From the 
electrostatic formulation of the Radon transform we 
can prove that the fidelity term in the functional Jms 
actually imposes that the electric field produced by the 
charge distribution / must be close to the “measured 
electric field”. Therefore we conclude that the term 
WTZf — g\\'^ 2 (p 3 -f is a fidelity term in this weaker sense. 

Using this property based on the electrostatic 
interpretation of the tomographic reconstruction, we 
can try to minimize some appropriate functionals in 
the new variables E (electric field) or ip (electric 
potential) and then compute the corresponding / 
(charge density). We manipulate the functional Jms 
as follows: 

>/Ms(r, /) 

= m - gWhiP^) +a [ |V/(a:)|2 dx + 

jR3\r 

= 2{2Tr)^\\E -Eg\\^ + a [ |V(V • E){x)\^ dx 
./R3\r 

+ /3'H^(r) 

= 2(27r)2||U -Eg\\^ + a [ \AE\^ dx + 

= EiT,E) (37) 

where U is a new functional depending on a vector 
function E and on a set T, and we used the fact that 
V A U = 0, since E is conservative. 

We observe that the functional 

E{T,E) = 2i2TTy\\E - EgW^ + a [ |AU|2 da; 

aH3\r 

+ I3H‘^{T) (38) 

is a second order functional for a vector field E 

in which appears the measure of the set F that is 

the set of discontinuities of / and thus is the set 
of discontinuities of V • U. In the functional E 
we recognize some similarities with a famous second- 
order free-discontinuity problem: the Blake-Zisserman 
model. This model is based on the minimization of the 
Blake-Zisserman functional 

BZ{ro,ri,v) = \\v-vo\\l2f_D) 

+ a |Ar;(a:)p dx 

./DXTouri) 

+ ^n^-^{TonD) 

+ jn^-\ir,\To)nD), ( 39 ) 


among admissible triplets (ro,ri,u), where 

• D C is an open set; 

• Fo,Fi c are closed sets; 

• Fq is the set of discontinuities of v (jump set), and 
Fi the the set of discontinuities of Vu (crease set); 

• v: D ^R,v G C^{D \ (Fq U Fi)) n C{D \ Fq) is 
a scalar function; 

• Av denotes the distributional Laplacian of v; 

• Vo & L'^{D) is the datum (grey intensity levels of 
the given image); 

• a, 13,^ > 0 are parameters; 

• denotes the {N — 1)-dimensional Hausdorff 
measure. 


The Blake-Zisserman functional allows a more precise 
segmentation than the Mumford-Shah functional in 
the sense that also the curvature of the edges of the 
original picture is approximated. On the other hand, 
minimizers may not always exist, depending on the 
values of the parameters /3 ,7 and on the summability 
assumption on vq . We refer to [29l [40l STJ |42l |43l |44] 
for motivation and analysis of variational approach to 
image segmentation and digital image processing. In 
particular see [301 ED SS] for existence of minimizer 
results and |42j for a counterexample to existence 
and oni on for results concerning the regularity of 
minimizers. 

Equation (37) implies that the functional Jms 
can be rewritten in terms of the vector field E and 
of the discontinuities set of V • U, i.e. the set of 
creases of E, using the terminology of the Blake- 
Zisserman model. The fact that in the functional F 
the discontinuities set of E is not present depends on 
the fact that we are assuming that the charge density 
/ in the functional Jms do not concentrate on surfaces 
or on lines. If we admit concentrated charge layers we 
can consider the Blake-Zisserman model for the vector 
function E as a relaxed version of the Mumford-Shah 
model for the charge /. In other words we propose to 
investigate the connections between minimizers of Jms 
and minimizers of the higher order functional Jbz- 

jBziro,ri,E) = \\E-Eg\\ L2(R3) 


+ a f |AU(x)p dx 

./R3\(rouri) 

+ /3H2(Fo) 

+ 7H2(ri\Fo), (40) 


with the additional constraint V A U = 0. The 
main advantage of this approach is that the functional 
Jbz is a purely differential functional, while the 
functional Jms is an integro-differential one. We 
expect that some results about the Blake-Zisserman 
model that could be rephrased into tomographic terms 
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would provide immediately new results in tomography. 
Conversely all the peculiar tomographic features as 
the intrinsic vector nature of the variable E, the fact 
that its support cannot be bounded and the extra¬ 
constraint V Ai? = 0, motivate new research directions 
in the study of free-discontinuities problems. For 
example, an interesting result in this context would 
be the determination of a good hypothesis on the 
datum Eg that ensure that the charge density / do 
not concentrate. 

We conclude this section with some comments: 


• We proved that the measured data g are actually 
the measured electric field produced by the 
unknown charge density, so the term \\TZf — 
511 ^ 2 (p 3 ) in the functional is a fidelity term in a 
weak sense. 

• The problem of the reconstruction of the charge 
can be rephrased into a reconstruction problem 
for the electric field. The electric field is an 
irrotational vector field, so the new minimization 
problem is actually a constrained minimization. 
In order to avoid this constraint one could 
reformulate the reconstruction problem in terms 
of the electric potential (p {E = — V:/?) obtaining 
a third-order functional in which the fidelity term 
is 


|jV(/7 — V(Pg|||2(R3), 

where the potentials are given by (36). 


(41) 


• All this considerations hold true in dimension 3. 
In a generic dimension n > 3 the situation 
is quite different because the inversion formula 
for the Radon transform involves a (possibly 
fractional) power of the Laplacian. In this case 
the electrostatic description of tomography given 
in this section fails. In order to restore it, 
it is necessary to consider another Radon-type 
transform which involves integrals of / over linear 
manifolds with codimension d such that (n — 
d)/2 = 1, i.e. d = n — 2, see e.g. [4} 1^. 
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